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Grant:  F49620-01-1-0390-P00003 

Title:  Numerical  Computation  in  MagnetoFluid-Dynamics 

Principal  Investigator:  Robert  W.  MacCormack 

Commercial  Phone:  650-723-4627  FAX:  650-725-3377 

Mailing  Address:  Department  of  Aeronautics  and  Astronautics 

_ Stanford  University 

_ Stanford.  CA  94305-4035 

E-Mail  Address:  nvmacc@aol.com 

AFOSR  Program  Manager:  Dr.  Fariba  Fahroo 

Research  Objectives: 

a)  Develop  algorithms  and  boundary  condition  procedures  for  solving  the  complete  equations 
governing  magneto-fluid-dynamics 

b)  Apply  the  developed  numerical  procedures  to  simulate  the  flow  about  realistic  aerospace 
configurations  of  Air  Force  interest 

c)  Support  and  complement  magneto-fluid-dynamic  research  studies  at  AFRL 
Status  of  Effort: 

Prior  to  the  start  of  the  subject  grant  the  following  research  accomplishments  were  made  by  the 
principal  investigator  toward  flow  simulation  governed  by  the  equations  of  magneto-fluid-dynamics. 

(1)  Reforming  the  flux  vectors  of  the  equations  governing  MFD  to  be  “homogeneous  of  degree 
one”  so  that  F  =  AU  and  a  flux  split  algorithm  could  de  devised.  F+  =  AJJ  and  F_  =  AU 

(2)  Simulation  of  the  flow  about  the  fore-body  of  a  hypersonic  vehicle  and  the  calculation  of  drag 
and  heat  transfer  with  and  without  magnetic  field  interaction. 

(3)  Splitting  the  flux  vectors  directly,  instead  of  basing  the  splitting  on  the  state  vector  U . 

F+  =  A+A~'F  and  F_  =  A_A~'F 

During  the  years  supported  by  the  subject  AFOSR  grant  the  following  research  accomplishments  were 
made. 

(4)  Inclusion  of  both  thermal  and  chemical  non-equilibrium  into  the  MFD  equations. 

(5)  Simulation  of  the  flow  within  MFD  generators  and  accelerators  for  the  proposed  energy 
“bypass  scram  jet  engine”  concept. 

(6)  Reformulating  the  governing  MFD  equations  for  strong  imposed  magnetic  fields,  of  the  order 
of  10  Tesla,  using  the  properties  of  the  imposed  magnetic  fields,  V-5  =  0  and  Vx 5  =  0 ,  to 
avoid  products  of  the  order  of  B 2  from  dwarfing  values  of  static  fluid  pressure. 

(7)  Analysis  of  the  physics  of  upstream  influence  caused  by  magnetic  diffusion  and  the  eigenvalue 
properties  of  the  MFD  equations  in  comparison  with  the  simpler  “low  magnetic  Reynolds 
Number”  approach  for  flows  within  MFD  generators  and  accelerators. 

(8)  Analysis  of  non-uniqueness  of  the  equations  of  magneto-fluid  dynamics 

(9)  Inclusion  of  Cesium  seeding  into  the  flow  and  the  Park  model  for  calculating  electrical 
conductivity  based  upon  electron  and  ion  concentrations. 

(10)  Analysis  of  the  physics  of  upstream  influence  of  the  flow  about  the  nose  of  a  hypersonic 

vehicle 


Accomplishments/New  Findings: 

Items  (4),  (6),  (8)  and  (9)  above  represent  progress  toward  research  objective  (a),  Items  (5)  and  (10) 
toward  research  objective  (b)  and  Items  (7)  and  (10)  toward  research  objective  (c).  Each  will  be 
discussed  below.  First,  we  present  the  governing  equations  of  magneto-fluid  dynamics. 

1.  The  Equations  of  Magneto-Fluid-Dynamics  - 

The  unsteady  equations  of  compressible  viscous  flow  within  an  imposed  magnetic  field  become 

dU  dF  dG  dH  . 

- +  —  +  —  +  —  =  0 

dt  dx  ay  dz 

The  state  flux  vector  is  given  by 

U  =  [p,pu,pv,pw,e  ,Bx,By,Bj\ 

with  density  p ,  velocities  u,v  and  w,  total  energy  per  unit  volume,  including  magnetic  field  energy, 
e’  =e  +  -^-B2 ,  and  Bx,  By  and  Bz  are  the  components  of  the  magnetic  field,  e  =  p(s  +  \{u2  +  v2  +  m’2))  , 

e  represents  the  internal  energy  and  B2  =  B2  +  B2  +  B2 .  The  flux  vector  F  becomes 
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with  p  =  p  +  -^-B2  and  magnetic  stress  given  by  j3v=-ve 


rdB±  dB^ 
dx,  dx. 


The  magnetic  field 


v  —t  ~"j  j 

components  shown  above  represent  the  total  of  the  imposed  and  induced  fields.  The  viscous  stress 


tensor  is  given  by  ttJ  =  -p 
vectors  G  and  H  are  similar. 


du,  | 

ydxj  dx,  J 


-8  A— where  8fJ  is  the  Kronecker  delta.  The  other  flux 
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TheB  JacobianB  ofB  theB  inviscidB  part®  ofB  theB  vectorB  F  B  (ie  .B  withB  p  —  A  —  Ve  —  0)B  isB 
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TheE  matrixE  AM  andE  eigenvaluesE  aboveE  wereE  givenE  byE  PowellE  alongE  withE  theE 
correspondingEeigenvectors .  H 

H 

TheE  f luxE  vectorsE  ofE  theE  MFDE  equationsE  areE  notE  homogeneous^  of B  degreeS  oneE  withE 
respectE  toE  theE  stateE  vector,  E  butE  theyE  canE  beE  triviallyE  modifiedE  soE  thatE  theyE 
becomeE  so.E  AE  methodB  forE  theirB  solutions  inE  conservations  lawE  formB  wasE  presentedB 


earlier3"4 .  E  TheE  methodE  usesE  upwindE  biasedE  f  luxE  vectorsE  directlyE  insteadE  of  E  usingE 
theB  stateE vectors . EThisB improvesE  accuracyE  andE  allowsE  f orE moreE f lexibilityE inB  theE 
eigenvectorEstructure , EwhichEwillEbeEexploitedElater . E 


Modified  Steger- Warming  flux  vector  splitting  is  used  to  approximate  the  inviscid  terms  and  central 
differences  are  used  for  the  viscous  terms.  A  block  implicit  algorithm,  using  modified  approximate 
factorization  with  sub-iteration,  was  then  used  to  solve  the  resulting  set  of  finite  volume  difference 
equations.  This  numerical  method  also  uses  frame  independent  smoothing,  and  has  the  option  of 
adapting  the  mesh  to  shock  waves.  The  gas  can  be  treated  as  a  perfect  gas  or  as  a  real  gas  in 
equilibrium. 

2.  Inclusion  of  Chemical  Non-Equilibrium  into  the  MFD  Equations  for  Simulation  of  the  Flow  within 
MFD  Generators  and  Accelerators  for  the  Proposed  Energy  “Bypass  Scram  Jet  Engine”  Concept. 
(Items  (4)  and  (5)  in  the  Status  of  Effort  Section  above). 

The  following  four  reaction  seeded  air  chemistry  model  of  Park,  Mehta  and  Bogdanoff,  for 
temperatures  to  4000K  was  used.  The  model  employed  here  contains  eight  species.  Species 
downstream  of  the  combustor  containing  hydrogen  were  not  included. 

1)  02  +  M  ^O+O  +  M 

2)  N2  +  0^>N0  +  N 

3)  0+N0<r*N+02 

4)  Cs  +  e~  <->  Cs+  +  e~  +  e~ 

The  forward  reaction  rates,  with  parameters  given  in  Tables  1-4,  are  of  fonn. 

kf(T)  =  AT”  exp(-7;  IT) 


Table  1.  02  +M  0  +  0  +  M 


me 

A,  0m8/  (kg-mol .  sec)  0 

E0N0 

Tr , EkE 

N2B 

BBEBE5 . 09xl0A12B 

-l.lE 

59, 360E 

02  E 

EEEBE5 . 09x10a12E 

-l.lB 

59, 360E 

NOE 

BBBEB1.27x10a13B 

-l.lE 

59, 360E 

NE 

EEEEE5 . 09x10a12B 

-l .  IE 

59 , 360B 

OE 

EBBBB5 . 09x10a12B 

-l.lE 

59 , 36 0E 

CsE 

BEEBES . 09x10 a12E 

-1 .  lE 

59 , 360B 

Cs+E 

EEBBES . 09x10 A12B 

-l.lE 

59 , 360B 

e-El 

E00005 . 09x10^140 

-1 .  lE 

59, 360E 

Table  2.  N2+0<r>N0+N 


A,  mV(kg-mol.sec) 

n 

Tr,  K 

5.69e6 

0.42 

42,938 

Table  3.  0+N0<^N  +  02 


A,  mVfkg-mol.sec) 

n 

Tr,  K 

2.36e3 

1.00 

19,220 

Table  4.C.S' +  e  f)Ci*+e  +e 


A,  ms/(kg-mol.sec) 

n 

Tr,  K 

3.90e27 

-3.78 

45,180 

The  backward  reaction  rates,  with  parameters  given  in  Table  5,  are  of  form 

kb(T)  =  kf(T)/keq(T)  , with  keq(T)  =  exp(AJ  z  + A,  + A,ln(z)  + A4z  + A,z2) 

and  z  =  10000/7\ 

Table  5.  Parameters  kea,  (MKS  units) 


Reac.  1 

Reac.  2 

Reac.  3 

Reac.  4 

A, 

-0.9278 

1.2441 

2.3074 

-4.1036 

a2 

17.1414 

0.7192 

-2.9933 

10.2278 

a3 

0.2816 

0.8606 

1.2493 

-3.4344 

a4 

-6.0607 

-3.9981 

-1.8594 

-4.2851 

A5 

0.0027 

0.00851 

0.0087 

-0.0097 

The  specie  mole  conversion  rate  equations  become 

Rt  ^{kf\02][Mm]-kbx[0}mMm}} 

m= 1 

R2=kf2[N2][0]-kb2[N0][N] 

R3=kf3[0][N0]-kb3[N][02] 

R4=kfA[Cs][e-]-kb4[Cs+][e-][e-] 

where  the  number  of  kilogram-moles/m3  of  species  Z  is  defined  by  [Z]  =  pj  o)2,  species  density  pz 
and  molecular  weight  co2 .  The  time  dependent  species  mole  equations  are  then 
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4 

The  heats  of  formation,  energy  taken  from  the  gas  internal  energy  to  form  species,  are  given  in  Table 
6.,  in  Joules/kilogram,  along  with  species  molecular  weights  and  the  initial  concentrations,  species 
density  divided  by  total  gas  density,  c2  =  pz!  p. 


Table  6. 1 

Physical  Properties  of  Gas 

c 0 

Initial  cz 

Heat  of  Form. 

N2 

28.01 

0.76290 

0.0 

02 

32.00 

0.23695 

0.0 

NO 

30.07 

0.0 

2.996123xl0A6 

N 

14.01 

0.0 

3.362161xl0A7 

o 

16.00 

0.0 

1.5431 19xlOA7 

Cs 

132.9 

0.00015 

0.0 

Cs+ 

132.9 

0.0 

2.899905xl0A6 

e- 

54.2* 

0.0 

0.0 

*  indicates  xlOA-5 


Also,  equations  for  species  convection  and  diffusion  need  to  be  solved. 

djuz  — 

dc.  dc ,  1  z  dx, 


dt  ‘  dxt  p  dx. 


z  =  1,8 


The  species  mole  equations  are  solved  implicitly  by  sub-iteration  in  time.  We  can  write  them  in  vector 

JO 

form  as  —  =  R(S ) .  S  is  a  vector  containing  the  8  species  mole  elements  and  R(S)  contains  the  8 
dt 


right  hand  sides  of  the  above  equations.  The  sub-iteration  procedure  converges  R(Sn+ )  - 


At 


to 


zero  as  follows. 

✓S  f  C(w)  _ 

{/  +  At  A}  (s(m+1)  -  S(m) )  =  At  R(S(m)) -  — 


with  Sm=S"  and  S(m)->Sn+l  as  At  is 


chosen  appropriately  small  and  A  =  — 1 — 

oS 


is  the  Jacobian  of  R  with  respect  to  S .  The  sub-iteration 


procedure  is  stopped  when  the  maximum  change  in  concentration  of  each  species  is  less  than  1 06 .  The 
species  convection-diffusion  equations  are  solved  using  upwind  implicit  finite  volume  approximations. 


2-D  simulations  of  the  flow  of  air  within  an  MHD  generator  were  made  with  seeded  cesium.  The 
generator  was  a  square  duct,  2.721m  long,  of  height/width  64.85cm  at  the  entrance  and  of  height/width 
85.14m  at  the  exit.  The  magnetic  field  across  the  channel  was  12.74  Tesla  and  the  transverse  voltage 
gradient  was  -29,400  V/m  at  the  entrance  and  -18,290  V/m  at  the  exit.  The  numerical  simulation  varied 
this  voltage  linearly  from  entrance  to  exit.  At  the  entrance  the  pressure  p=l,038xl06  Pa,  temperature 
T=3371°K,  and  the  Mach  number  equaled  2.418.  The  electrical  conductivity  was  enhanced  by  the 
seeding  with  liquid  cesium.  The  Park  program  for  calculating  conductivity  had  difficulties  near  the 
wall,  which  was  held  at  a  fixed  wall  temperature  of  T=300°K,  and  had  to  numerically  limited  within 
the  boundary  layer  to  levels  predicted  by  the  Park  program  at  the  edge  of  the  boundary  layer.  This 
difficulty  was  later  removed,  as  will  be  discussed  later.  The  results  for  finite  rate  chemistry  and 
conductivity  calculated  from  the  Park  program  are  shown  in  Figs.4-6. 


B 

Figure  1.  Velocity  profiles  within  MHD  generator,  reacting  seeded  air  chemistry.B 


TheBsupersonicBvelocityBprof ilesBwithinBtheBgenerator , BshownBalongBtheBcenterB 
planeBinBFig . 1 , BdecelerateBsignif icantlyBinBtheBdivergingBductBwhileBgeneratingB 
electricBpower . B 


Figure  2.  Normalized  force  and  velocity  in  axial  direction  within  MHD  generator  section,  reacting 
seeded  air  chemistry 

NormalizedBaxialBforceBFjBandEmaximunnBvelocityBasBaBfunctionBof  B  x  BareBshownBinE 
Fig .  2  .  ENoteBthatBwithinEtheBgeneratorBbothB  Z^BandB  «max  BdecreaseBsignif  icantly .  B 


3.  An  Alternative  Formulation  of  the  Equations  of  Magneto-fluid-Dynamics  -  The  Reduced  MFD 
Equations  in  Conservation  Form  (Item  (61  in  the  Status  of  Effort  Section  above). 

An  alternate  formulation  of  the  governing  equations  of  magneto-fluid-dynamics  has  been  devised.  This 
formulation  is  mathematically  equivalent  to  the  original  set  of  equations  governing  magneto-fluid- 
dynamics  given  above,  including  magnetic  self  induction,  and  retains  the  conservation  law  form  of  the 
equations  and  their  eigenvalues.  This  new  set  has  advantages  over  the  original  set  when  solved 
numerically  for  flows  within  strong  imposed  magnetic  fields. 

The  formulation  treats  the  imposed  and  induced  fields  separately.  Though  ever  present,  the  imposed 
field  remains  in  the  background  as  the  finite  volume  difference  equations  focus  on  the  induced  field. 
The  imposed  field  can  not  be  eliminated  entirely  from  the  difference  equations  because  of  non¬ 
linearity,  but  no  squared  imposed  magnetic  field  terms  appear.  This  is  important  for  flows  within 
strong  imposed  magnetic  fields,  because  the  magnetic  pressure,  proportional  to  the  square  of  the 
imposed  magnetic  field  strength,  can  be  several  orders  of  magnitude  larger  than  the  aerodynamic 
pressure  or  the  induced  magnetic  field  pressure.  Numerical  errors  in  the  very  large  magnetic  stress 


difference  terms  of  the  imposed  field  could  be  of  significance  when  combined  with  the  relatively 
smaller  fluid  stress  terms. 

Before  the  introduction  of  this  alternative  formulation,  there  were  two  choices  for  including  the  effects 
of  a  magnetic  field  upon  an  ionized  flow:  (1)  the  complete  equations  of  magneto-fluid-dynamics, 
including  magnetic  induction  and  (2)  the  inclusion  of  the  “j  cross  B”  force  and  Joule  heating  effects  in 
the  Navier-Stokes  equations  as  additional  source  terms.  The  first  choice  is  a  set  of  eight  equations 
consisting  of  the  Navier-Stokes  equations,  with  added  magnetic  stress  tensor,  plus  the  Maxwell 
inductions  equations  and  is  presented  in  Sec.l.  above.  The  second  choice,  a  set  of  five  equations, 
called  the  ’’Low  Magnetic  Reynolds  Number  Approximation”,  assumes  that  the  induced  magnetic 
field  is  negligible.  It  is  far  more  efficient  and  has  far  fewer  numerical  difficulties  associated  with 
inclusion  of  magnetic  effects  than  the  first  choice.  There  is,  however,  some  uncertainty  in  when  the 
low  magnetic  Reynolds  number  is  valid,  even  for  aerodynamic  flows  of  current  interest.  The  now 
available  third  choice  presented  below  can  be  used  to  investigate  this  question. 

The  total  magnetic  field  consists  of  the  imposed  magnetic  field,  Bo ,  and  the  induced  magnetic  field, 
Bi,  and  Bt  =  Bo  +  B, ,  where  the  subscripts  t,  0  and  /  now  and  below  indicate  total,  imposed  and 
induced  magnetic  components.  For  cases  for  which  the  induced  field  is  much  less  than  the  imposed 

field,  but  not  negligibly  small,  we  can  benefit  by  rewriting  the  Lorentz  force  as  Z/  = — (VxB,)xft , 

Me 

because  the  imposed  magnetic  field  is  generated  by  currents  external  to  the  flow  field,  for  which 
Vx5o  =  0.  The  approach  taken  here  is  similar  to  the  simplification  in  electromagnetic  scattering 
where  only  the  disturbed  field  is  calculated,  with  nothing  lost  by  the  separation  of  the  two  fields. 

We  can  also  write  the  Lorentz  force  as  Z/=  —  (Vxfi)xA — -(VxiFo)xZo,  and  through  some 

Me  Me 

algebraic  manipulation,  the  Lorentz  force  can  be  brought  into  the  flux  derivative  terms  of  the 
momentum  equation,  in  conservation  law  form  as  before.  The  state  vector  becomes 

U  =  [a  pu, pv, pw, e' ,  B,x ,  Biy,  Biz ]  and  the  new  flux  vector  F  becomes 


F  = 


pu 

pu2  +  p  +  --^BixBixyv--^{BixBox} 
pvu  +  rxy  -jrBiyBtx  -±[BixBoy\ 
pwu  +  Txz  -jrt  BizBt  -  f'  {BixBoz } 

0  J1 

(e  +P+TJU  +  txyv  +  TXZW  -k— 

+j-(-(Bi  ■  u)Bix  +  PxxBix  +  PxyB,y  +  &A ) 


with 


A, 

B, yu-B,xv  + pxy 
B,zu-B,xw  +  fixz 


dBij 

dx. 


dx 


J  J 


J 


p  =  p  +  jpB,2  +-±-Bi  Bt  and  the  magnetic 


stress  given  by  j3tJ  =  -ve 


By  replacing  the  magnetic  pressure,  ^B,2,  by  the  smaller  -±-B2  +±B,B,,  the  magnetic  and  static 

pressures  are  closer  in  magnitude  for  strong  imposed  magnetic  fields.  Favorable  reductions  also  take 
place  in  the  induction  equations  because  the  magnetic  diffusion  terms  j3;J  are  just  components  of  the 

curl  of  the  induced  magnetic  field  times  ve .  Again  the  imposed  field,  produced  by  currents  outside  the 
flow  field,  is  curl  free.  Hence,  veV xB,  =  veV xBi .  Here  also  the  production  and  diffusion  terms  are 
more  equally  balanced.  Finally,  the  terms  in  the  curly  brackets  in  the  momentum  equations  above 
vanish  if  the  imposed  field  is  constant  in  space  because  of  the  divergence  free  nature  of  both  the 
imposed  and  induced  fields. 

One  may  assume  that  the  structural  changes  to  the  equations,  just  presented,  from  the  separation  of  the 
induced  and  imposed  fields,  would  have  profound  changes  to  the  original  eigenvalue-eigenvector 
structure  of  the  equations.  Fortunately,  the  eigenvalues  remain  the  same  as  shown  below. 
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The  eigenvectors  are  changed,  however,  but  the  original  set  can  still  be  used  in  the  solution  procedure 
as  is  to  solve  the  alternative  RMFD  (Reduced  Magneto-Fluid  Dynamics)  equations  just  presented, 
because  of  the  conservation  form  of  the  flux  vector  splitting  used. 


The  term  Reduced  is  used  here  to  reflect  the  notion  that  the  magnitude  of  the  magnetic  terms  are 
reduced  by  removing  as  often  as  possible  the  imposed  magnetic  field  from  them,  although  the  number 
of  terms  is  actually  increased.  The  RMFD  equations  are  mathematically  and  physically  equivalent  to 
the  original  MFD  equations. 


4.  The  Navier-Stokes  Equations  Plus  a  Source  Term  for  Magnetic  Field  Effects  -  The  Low  Magnetic 
Reynolds  Number  Approximation 

An  ionized  flow  within  an  imposed  magnetic  field  can  self  induce,  thus  changing  the  magnitude  of  the 
total  magnetic  field.  The  relative  magnitude  of  the  induced  component  depends  upon  the  Magnetic 
Reynolds  number,  defined  by  Rm  =  u0l0crejue ,  where  w0  and  /0  are  reference  flow  speed  and  length,  <7e 

is  the  gas  conductivity  and  ,  fxe  is  the  magnetic  permeability.  The  magnetic  diffusion  coefficient  is 

given  by  ve  =  — .  For  most  aerodynamic  flows  the  gas  conductivity  is  very  small  and,  consequently, 

ve  is  very  large  and  the  magnetic  Reynolds  number  is  less  than  one.  In  such  cases,  any  self  induced 
magnetic  field  supposedly  rapidly  diffuses  away,  leaving  only  the  imposed  magnetic  field.  This  leads 
to  a  great  simplification  in  calculating  the  electro-magnetic  effects  upon  an  ionized  flow.  This 
approach,  The  Low  Magnetic  Reynolds  Approximation  approach ,  only  needs  to  add  a  source  term  to 
the  flow  equations  given  above  to  include  the  electro-magnetic  field  effects. 


Ionized  flow  in  the  presence  of  a  magnetic  or  electric  field  generates  a  current,  according  to  Ohm's 
law,  j  =  ae(E  +  uxB),  where  j  is  the  current  density,  E  is  the  electric  field  potential,  u  is  the  flow 

velocity  and  B  is  the  magnetic  field.  The  electric  current  itself  interacts  with  the  magnetic  field  to 
create  a  Lorentz  force,  L /  =  jxB ,  that  acts  on  the  flow  in  addition  to  pressure,  p,  and  viscous  stress. 


In  addition  to  the  Lorentz  force  added  to  the  momentum  equations,  Joule  heating,  caused  by  the  flow 
of  electric  current  through  the  fluid,  plus  magnetic  force  work  terms  need  to  be  added  to  the  energy 
equation.  The  equations  become 

dU  dF_  dG  dH_  =  s 
dt  dx  dy  dz 

The  state,  flux  vectors  and  the  source  vector  are  given  by  U  =  \p,pu,pv,pw,e\ 


( 


F  = 


pu 

pu2  +  p+r„ 
pvu  +  rxy 

pWU  +  Tx: 

0 

(e  +  p  +  Tju  +  Txyv  +  T^W  -k—j 


,  etc,  HHandHHH5  = 


0 

(7  xB)x 
(7  xB)y 

(jxB)z 

0  xB)-u  +  -^J  ■] 


This  equation  set  is  sufficient  to  describe  ionized  flow  within  an  electro-magnetic  field,  as  long  as  the 
fields  are  specified.  It  is  not  much  more  difficult  to  solve  than  the  underlying  Navier-Stokes  flow 
equations  themselves.  However,  if  the  magnetic  field  varies  in  time  by  self  induction  and  the  induced 
magnetic  components  are  relatively  significant  in  magnitude  then  the  equations  for  magnetic  induction 
also  need  to  be  solved.  This  larger  set,  shown  earlier  in  Sec.  1 .,  is  much  more  difficult  to  solve  and 
should  not  be  attempted  if  it  can  be  avoided. 


5.  Comparison  of  the  Full  MFD  Equations,  the  Reduced  MFD  Equations  ans  the  Low  Reynolds 
Number  Equations 

The  figure  below  shows  the  velocity  field  within  the  accelerator  region  of  a  proposed  energy  bypass 
scram  engine,  locate  just  downstream  of  the  combustor. 


Figure  3.  Velocity  profiles  within  an  MHD  accelerator,  reduced  MHD  simulation 


E 

Figure  4.  Normalized  Force  and  velocity  in  the  axial  direction  within  an  MHD  accelerator  section,  (a) 
Full  MFD,  (b)  Reduced  MFD  and  (c)  Low  Magnetic  Reynolds  Number 

The  supersonic  velocity  profiles  within  the  accelerator  are  shown  in  Fig.3  to  accelerate  significantly  in 
the  converging  duct.  These  results  were  obtained  from  the  Reduced  MFD  set  of  equations.  This  set  of 
equations  is  mathematically  equivalent  to  the  full  set  of  MFD  equations,  but  the  imposed  and  induced 
magnetic  fields  are  separated.  The  imposed  field  is  constant  and  can  often  be  extracted  from 
derivatives  in  the  governing  equations.  This  is  numerically  helpful  when  the  imposed  field  is  much 
larger  than  the  induced  magnetic  field  and  consequently  the  imposed  magnetic  pressure,  proportional 
to  the  magnitude  of  the  imposed  field  squared,  is  orders  of  magnitude  greater  than  the  aerodynamic 
pressure.  The  reduced  set  of  equations,  called  Reduced  MFD  herein,  is  still  in  conservation  form  and 
has  the  same  eigenvalue  structure,  but  is  numerically  better  posed.  The  Low  Magnetic  Reynolds 
Number  Approximation  neglects  the  induced  terms  entirely.  The  Reduced  MFD  formulation  was  not 
foreseen  in  the  original  proposed  research,  but  developed  subsequently  after  strong  magnetic  fields 
were  introduced  into  the  study  to  simulate  flows  of  Air  Force  interest. 

Normalized  axial  force  Fx  and  maximum  velocity  wmax  are  shown  within  the  accelerator  in  Fig.4.  All 

three  simulations,  (a)  the  full  MFD,  (b)  the  reduced  MFD  and  (c)  the  Low  Magnetic  Reynolds  Number 
Approximations  are  shown.  The  induced  magnetic  field  for  both  the  Reduced  and  full  MFD 
approaches  was  approximately  1%  of  the  imposed  field  strength.  Note  for  each  case  that  both  Fx  and 

Mmax  increase  significantly  at  first  then  level  off  as  the  flow  traverses  the  duct.  The  full  MFD  results  for 
velocity  show  unexpected  spatial  variations,  perhaps  caused  by  the  numerical  difficulties  discussed 
earlier.  Also,  the  Low  Magnetic  Reynolds  Number  Approximation  is  not  in  agreement  for  calculated 
thrust  with  the  other  results  that  include  magnetic  induction.  It  is  about  20%  lower.  This  finding  can  be 
very  significant.  If  the  Low  Reynolds  Number  Approach  is  found  to  be  physically  incorrect  for  flows 
at  magnetic  Reynolds  Numbers  less  than  one,  as  in  the  case  here,  their  use  to  simulate  flows  within 
magnetic  and  electric  fields  will  have  little  confidence.  The  small  1%  change  in  magnetic  field  caused 


by  induction  results  in  a  magnetic  pressure  of  the  order  of  the  gas  pressure  itself  and  could  explain  this 
effect.  Further  study,  required  to  determine  the  significance  and  validity  of  these  results  is  pursued  in 
the  next  section. 


6.  Analysis  of  the  Physics  of  Upstream  Influence  Caused  by  Magnetic  Diffusion  and  the  Eigenvalue 
Properties  of  the  MFD  Equations  in  Comparison  with  the  simpler  “Low  Magnetic  Reynolds  Number” 
Approach  for  Flows  within  MFD  Generators  and  Accelerators.  (Item  (7)  in  the  Status  of  Effort  Section 
above). 

6. 1  Simple  test  geometry  problem 

TheseS  figures®  demonstrated  that®  the®  two®  mathematical®  formulations®  described 
dif  ferentd  physics  .  d  Ond  thed  Shownd  below®  ared  thed  resultsd  f  romd  thed  threed  numericald 

simulations® for® flow® within® an® MFD® accelerator® with® a® simplified® geometry . ® The® 2D® 

channeld  consistsd  ofd  ad  rectangulard  volumed  3m. d  longd  andd  lm.d  high.d  Thed  walld 

temperaturesdweredhelddf  ixeddatd3  00°K.  dThedimposeddmagneticdf  ielddwasd  By  =  lldTeslad 

andd  thed  electricd  fieldd  wasd  Ez  =-30,  OOOd  V/m.d  Thed  eguilibriumd  flowd  enteringd  thed 
channeld  wasd  atd  Machd  1 . 166d  withd  totald  pressure ,  d  8xl05d  N/m2d  andd  totald  temperatured 
7500°K.d  Thed  electricald  conductivityd  wasd  heldd  uniformd  withind  thed  acceleratord  atd 

O  =  36.0 /(Sim)  ®and®the®magnetic®Reynolds®number® =  U0l0CtePe  =  0.1847  . 
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The  acceleration  of  the  velocity  vectors  is  shown  in  Fig.5.  These  results  were  obtained  from  the 
Reduced  MFD  set  of  equations.  Normalized  axial  force,  Fx(x )  =  J  (pu2  +  p)  ds ,  integrated  across  the 

S(x) 

channel,  and  maximum  velocity  wmax  are  shown  within  the  accelerator  in  Fig.6  for  all  three 
simulations.  The  maximum  induced  magnetic  field  was  less  than  2%  of  the  imposed  field  strength. 
Note  for  each  case  that  both  Fx  and  wmax  increase  significantly  at  first  then  level  off  as  the  flow 
traverses  the  duct.  The  Full  MFD  results  for  velocity  show  some  unexpected  spatial  variations,  perhaps 
caused  by  the  numerical  difficulties  discussed  earlier.  Also,  the  Lowi?m  MFD  result  is  not  in 

agreement  for  calculated  thrust  with  the  other  results  that  include  magnetic  induction.  It  is  about  8% 
lower.  The  small  2%  change  in  magnetic  field  caused  by  induction  results  in  a  magnetic  pressure  of  the 
order  of  the  gas  pressure  itself  and  was  first  thought  to  be  responsible  for  this  effect 


Figure  5.  Velocity  profiles  within  an  MHD  accelerator,  Reduced  MHD  simulation 


Figure  6.  Normalized  Force  (bottom)  and  velocity  (top)  in  the  axial  direction  within  an  MHD 
accelerator  section,  (a)  Full  MFD  -  diamonds,  (b)  Reduced  MFD  -  circles  and  (c)  Low  Rm  MFD  - 

squares 

6.2  Possible  Causes  for  Differences 

Consider  the  following  explanations  for  the  differences  in  computed  results: 

(1)  Errors  in  the  numerical  procedures  caused  the  difference  in  results. 

(2)  The  small  induced  magnetic  fields  are  significant  enough  to  change  the  magnetic  stress,  work 
and  heating  within  the  flow  field. 

(3)  The  eigenvalue-vector  structure  of  the  equations  is  different  enough  to  significantly  alter 
domains  of  dependence  and  influence  within  the  flow  field. 

6.2. 1  Numerical  error  in  solving  the  equations 

Unfortunately,  numerical  error  is  always  present  despite  the  most  exhausting  steps  to  prevent  it  and  it 
should  be  the  first  consideration.  The  test  channel  accelerator  problem  was  made  as  simple  as  possible 
to  reduce  potential  error.  The  geometry  is  Cartesian-like,  the  magnitudes  of  the  magnetic  field,  electric 
field  and  electrical  conductivity  were  constants,  the  flow  was  supersonic,  except  within  the  boundary 
layer  and  the  solutions  were  steady  state.  A  single  program  was  written  to  solve  each  set  of  equations. 


Hopefully,  only  the  equations  changed  and  not  the  solution  procedure  for  solving  them.  The  Low  Rm 

MFD  set  of  equations  was  always  the  easiest  to  solve  and  the  fastest  to  converge.  The  Full  MFD 
equations  were  the  most  difficult  to  solve,  but  they  should  give  the  right  answer  if  solved  correctly  and 
should  be  able  to  be  used  to  tell  when  the  simpler  sets  suffice.  However,  it  will  probably  be  quite 
awhile  before  we  can  solve  them  with  confidence  for  all  engineering  flows  of  interest  to  hypersonic 
vehicle  design.  Much  more  experimentation  will  be  required  for  validation  of  numerical  simulations  of 
MFD  flows. 


Figure  7.  Normalized  Force  and  velocity  (a)  Reduced  MFD  -  circles,  (b)  original  Low  R m  MFD,  - 
squares  and  (c)  Low  Rm  MFD,  using  the  total  magnetic  field  from  the  above  Reduced  MFD  solution 
-  small  squares 

6.2.2  Induced  magnetic  field  may  significantly  change  magnetic  stress 

To  test  out  this  hypothesis  an  additional  calculation  was  made.  The  induced  magnetic  field  from  the 
solution  of  the  Reduced  MFD  set  of  equations  was  added  to  the  imposed  field  and  used  as  an  initial 
condition  for  the  Low  Rm  MFD  set  of  equations.  Therefore,  both  sets  of  equations  had  the  same 
magnetic  fields  at  steady  state.  The  results  are  shown  in  Fig.7.  Although  there  are  small  differences 
between  the  two  Low  7^  MFD  solutions  they  end  up  with  the  same  thrust  and  velocity  at  the  exit. 
Therefore  the  small  difference  in  magnetic  induction  is  not  the  cause  of  the  difference  in  the  thrust 
observed  earlier. 

6. 2. 3  The  different  eigenvalue-vector  structure 

The  eigenvalue-vector  structure  of  the  equations  is  different.  The  Low  Rm  MFD  set  of  five  equations 
has  the  same  eigenvalues  and  eigenvectors  of  the  underlying  Euler  or  Navier-Stokes  equations.  These 
eigenvalues  for  the  flux  in  the  x-direction  are  \  34  =  u  and  A25  =  u±c,  where  u  is  the  velocity  in  the 

x-direction  and  c  is  the  speed  of  sound.  These  eigenvalues  determine  the  domains  of  mathematical  and 
hence  assumed  physical  domains  of  dependence  and  influence.  For  the  MFD  channel  accelerator  case 
the  flow  is  supersonic  except  within  a  very  thin  boundary  layer  at  the  walls.  Therefore  all  the 


eigenvalues  are  positive  and  permit  no  upstream  influence.  The  Navier-Stokes  equations  also  permit 
diffusion  because  of  viscosity.  But  outside  the  boundary  layer  viscous  diffusion  is  insignificant  for  the 
flow  Reynolds  number  of  the  test  problem. 

The  eigenvalues  for  the  flux  in  the  x-direction  for  the  Full  or  Reduced  MFD  eight  equation  set  are 
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where  B2  is  the  square  of  the  magnetic  field,  Bx  is  the  x-component  and  p  is  the  density.  are  the 
Alfven  wave  speeds  and  and  AyS  are  called  the  fast  and  slow  magnetic  wave  speeds.  If  the 
magnetic  field  is  zero  in  magnitude,  A,  8  vanish  and  both  sets  of  eigenvalues  are  equal.  In  addition  to 

viscous  diffusion,  magnetic  diffusion,  which  can  be  orders  of  magnitude  larger,  can  spread  information 
omni-directionally.  Even  though  the  flow  is  conventionally  said  to  be  supersonic  in  the  x-direction  for 
the  test  case  problem,  all  the  eigenvalues  are  not  positive  and  information  can  easily  travel  upstream. 
The  magnetic  Mach  number  Mmag  =  u/cf=  0.148  at  the  entrance,  where  cf  is  the  fast  magnetic  wave 

speed  and  is  given  above  by  the  large  square  root  term  appearing  in  the  equation  for  A^  above. 

The  governing  equations  are  meant  to  match  the  assumed  physics  and  moreover  vice  versa.  Here  the 
physics  corresponding  to  the  two  descriptions,  Full  MFD  and  Low  Rm  MFD  are  definitely  different.  Is 
this  difference  responsible  for  observed  differences  in  the  flow  simulations? 

The  answer  is  yes,  which  we  can  demonstrate  with  a  little  analysis.  First,  the  Lorentz  force  is  given 
below. 


6.3  The  different  eigenvalue-vector  structure 

The  eigenvalue-vector  structure  of  the  equations  is  different.  The  Low  Rm  MFD  set  of  five  equations 


Lf  =—(VxB)xB  =  o-  (E  +  uxB)xB .  The  LowR„,  MFD  equations  evaluate  the  Lorentz  force  term 
Me 

by  the  expression  to  the  right  of  the  last  equal  sign  above  and  the  Full  MFD  equations  by  that  to  the 
right  of  the  first  equal  sign  above.  For  the  two  dimensional  test  flow  problem  the  curl  of  the  magnetic 
field  has  only  a  z-direction  component. 


A  =  (Vx5,)z  = 


dy 


Second,  we  can  express  the  induction  equations  as 

dBx  _  d(vBx  -  uBy  +  veDz )  _  BE, 

dt  dy  dy 


dBy  |  d(vBx-uBy  +  veDz)  _  ,  dEz 

dt  dx  dx 

and  finally  from  them  form  a  Poisson  equation  for  Dz  as  follows. 

dTf  _  d%  d2Ez 

dt  dx2  dy2 


We  can  solve  this  equation  in  an  extended  domain  of  the  channel  flow  problem  to  illustrate  the 
powerful  upstream  and  downstream  influence  of  the  induction  equations.  Earlier  the  domain  enclosed 
only  the  flow  between  the  two  electrodes  extending  from  x=0  to  x=3m.  The  extended  domain  covers 
the  flow  ahead  and  downstream  of  the  electrodes  by  lm.  each. 

Figs.8  and  9  show  the  solutions  for  Ez  and  D:  from  the  Poisson  equation  on  this  domain.  They 
represent  solutions  to  the  Full  and  Reduced  MFD  equations. 


Figure  8.  Contours  of  Ez  on  extended  domain 


Figure  9.  Contours  of  Dz  on  extended  domain 

Note  the  upstream  propagation  of  Ez  for  the  Full  MFD  solution  in  Fig.8  and  the  strong  gradients  and 
upstream  propagation  in  Dz  shown  in  Fig. 9,  both  of  which  could  not  be  simulated  using  the  Low  Rm 
MFD  approach. 

6.4  Conclusion  for  test  problem  results 

These  figures  demonstrate  that  the  two  mathematical  formulations  describe  different  physics.  On  the 
original  domain  the  Low  Rm  MFD  description  behaved  exactly  as  a  conventional  supersonic  flow  with 

no  apparent  upstream  influence.  On  the  other  hand,  the  Full  MFD  approach,  through  the  fast  magnetic 
wave  speed  and  high  magnetic  diffusion,  did  allow  upstream  influence  even  at  the  upstream  boundary, 
from  which  the  flux  split  procedure  accepted  only  information  carried  along  characteristics  needed  for 
the  downstream  solution.  The  differences  shown  earlier  in  thrust  prediction  resulted  not  from  different 
velocity  distributions  but  from  different  mass  flow  rates  entering  the  channel  caused  by  upstream 
influence. 

However,  there  is  considerable  doubt  that  placing  this  boundary  at  the  start  of  the  electrodes  for  the 
Full  MFD  set  of  equations,  as  shown  above  for  the  simple  test  case  geometry,  allowed  proper  upstream 
influence  in  view  of  the  rapid  gradients  shown  above  at  x— Om.  in  Fig.9.  These  calculations  need  to  be 
repeated  with  the  upstream  boundary  located  sufficiently  far  ahead  of  the  electrodes,  as  shown  in  the 
next  subsection. 

6.5  Extended  accelerator  geometry  test  problem 

The  channel  was  extended  1 .5  meters  both  ahead  and  aft  of  the  imposed  electrode  and  magnetic  field 
section,  thereby  placing  the  boundaries  at  sections  of  the  flow  that  are  supersonic  with  respect  to  all 
wave  speeds  of  sound.  The  equations  of  the  Full  MFD  Approach  will  induce  a  magnetic  field  beyond 
the  electrode-imposed  magnetic  field  location,  but  it  should  not  be  strong  enough  to  produce  magnetic 


waves  capable  of  reaching  the  upstream  boundary.  The  flow  within  the  extended  channel  will  be 
simulated  again  using  both  approaches  and  then  analyzed  and  discussed. 
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Figure  10.  Sketch  of  Extended  MFD  Accelerator 


Figure  11.  Imposed  Magnetic  Field 


Figure  12.  Imposed  Electric  Field 


The  extended  MFD  accelerator  geometry  is  shown  in  Fig.  10  and  the  imposed  electric  and  magnetic 
fields  in  Figs.  1 1  and  12.  The  discontinuities  in  these  fields  present  simulation  difficulties.  For  example, 
the  magnetic  pressure  suddenly  jumps  to  a  value  several  times  that  of  the  fluid  static  pressure  at  x=0 
and  then  back  down  again  at  x=3m.  The  equations  describing  this  flow  and  a  needed  reformulation  for 
them,  designed  to  overcome  the  simulation  difficulty,  are  given  in  the  following  section. 


The  imposed  magnetic  field  was  Byo  =4  Tesla  and  the  electric  field  was  Ezo  --10,909  V/m.,  creating  a 

load  factor  of  2.004.  The  wall  temperatures  were  held  fixed  at  300°K  and  the  initial  flow  speed  was 
1361m/sec  The  equilibrium  flow  entering  the  channel  was  at  Mach  1.166,  with  pressure  1.251x10 
N/m2  and  temperature  3,583°K.  The  electrical  conductivity  was  held  uniform  within  the  accelerator  at 
a  =  36.0 /(Qw)  and  the  magnetic  Reynolds  number  Rm  =  u0l0aeMe  =  0. 1 847 . 


.The  solutions  for  By  and  Ez  are  shown  in  Figs.  13  and  14,  along  the  channel  centerline.  Note  the 
relatively  small  induced  magnetic  field.  The  solution  of  the  ’’Equation  for  Ez  ”  of  Subsec.6.3  above  is 
denoted  by  E\ .  It  was  solved  independently  of  the  MFD  equations  of  Subsec.6.3.  The  solution  of  the 
MFD  equations  is  denoted  by  Ez  and  the  solid  symbols.  The  two  solutions  are  in  excellent  agreement, 
which  is  a  validation  of  the  MFD  simulation. 


Figure  15.  Bx  (top),  By  (middle)  and  Ez  fields. 
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Figure  16.  Velocity  Field 
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Figure  17.  Thrust,  Fx(x )  =  J  ( pu 2  +  p)  ds 

S(x) 


6.5.1  A  Serious  Problem  Concerning  the  Results  Given  Above 

The  solution  to  the  MFD  Equations  was  fully  converged.  That  is,  it  satisfies  the  governing  equations 
relaxed  to  steady  state  and  it  also  is  in  excellent  agreement  with  the  results  for  the  ’’Equation  for  Ez  ” 
of  Subsec. 6.3  above.  But,  perhaps  the  reader  noticed  that  it  does  not  satisfy  the  constraint  that  the 
induced  field  be  divergence  free. 

6.6  Divergence  ofB  Control 

6. 6. 1  Artificial  “Compressibility  ”  Approach 

For  incompressible  flow,  artificial  compressibility  is  often  used  to  steer  the  velocity  field  toward  being 
divergence  free.  An  artificial  compressibility  term  proportional  to  the  divergence  of  the  velocity  field 
is  added  to  the  pressure,  p<r-  p  +  ■  u) .  Similarly,  to  control  divergence  of  the  magnetic  field, 

we  add  the  term  ve(V-B)  to  the  induction  equations. 

New  2-D  Induction  Equations  - 

dBx=  dve(VB)  d{vBx-uBy+veDz) 
dt  dx  3 y 

dB}  |  d(vBx  - uBy  +veDz)  d (V •  B) 
dt  dx  dy 

The  results  given  earlier  did  use  the  above  equation,  but  it  was  not  sufficient  to  prevent  V  •  B  from 
becoming  too  large.  Though  it  did,  perhaps,  prevent  runaway  growth. 


Partial  Poisson  Equation  Approach 


d2<j) 

A  1-D  Poisson  Equation  for  <p  is  solved  at  each  flow  field  point. 

component  is  kept  as  is,  but  the  y  component  is  then  modified, 
divergence  free  magnetic  field. 


=  -V  •  B .  The  x  magnetic  field 
dd> 

B  <r- B  +  — ,  to  maintain  a 
'  '  dy 


6. 7  New  Computational  Results  for  Extended  MFD  Accelerator 

The  flow  field  was  recalculated,  using  the  Partial  Poisson  Equation  Approach  discussed  above  and  the 
new  results  are  shown  in  the  following  figures. 


Figure  18.  Solution  for  By  along  centerline. 

Notice  that  By  now  increases  at  x=0,  where  previously  it  decreased,  as  shown  in  Fig.  13. 


Figure  19.  Solutions  for  Ez  along  centerline. 

Fig.  19  reveals  a  really  shocking  surprise.  The  Ez  solution  to  the  MFD  equations,  with  the  constraint  on 
the  divergence  of  the  magnetic  field  maintained,  adjusts  toward  the  original  imposed  magnetic  field 
EZOj  and  not  the  solution  E\  to  the  ’’Equation  for  E,  ”  of  Subsec.6.3  given  earlier. 


Without V -B  =  0  control  by  a  Poisson  Eq. 
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With  V  •  5  =  0  control  by  Partial  Poisson  Eq. 
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Figure  20.  Magnetic  Field  Lines 


The  magnetic  field  lines  for  both  sets  of  results  are  shown  in  Fig.20.  The  results  of  Sec.6.5  did  not  use 
a  Poisson  Equation  to  maintain  a  divergence  free  magnetic  field  and  the  results  of  the  present  section 
used  a  Partial  Poisson  Equation  for  control.  Note  the  large  difference  in  the  field  line  distributions.  The 
top  figure  resembles  a  supersonic  wave  pattern  in  the  electrode  region  while  the  lower  figure  appears 
to  be  like  streamlines  for  a  slowly  rotating  flow. 


Effort  Section  above). 

Two  solutions  have  been  presented  for  the  extended  accelerator  geometry  test  problem.  Is  either 
solution  correct?  Both  satisfy  the  MFD  equations  and  boundary  conditions.  Compare  the  solutions  for 


Ez  in  Fig.  14  with  that  in  Fig.  1 9.  To  say  that  the  latter  solution  is  correct  because  it  satisfies  V  •  B  =  0  is 
not  sufficient.  The  first  solution  agrees  with  the  solution  E\  to  the  ’’Equation  for  E”  of  Subsec.6.3 


and  the  latter  adjusts  toward  the  original  imposed  magnetic  field  Ew  •  We  can  determine  which  is 
correct  by  looking  at  the  unsteady  Ampere-Maxwell  Equation. 


BE_  1  T  VxB  j 

dt~  e\  Me  J 

But  an  underlying  assumption  made  in  deriving  the  MFD  Equations  is  J  —  V  X  B 

Hence  —  =  0  and  E  =E  is  the  correct  solution.  It  satisfies  the  unsteady  Ampere-Maxwell 
’at  z  ” 

Equation  as  well  as  the  constraint  V  •  B  =  0 . 


The  velocity  and  By  fields  are  shown  in  Fig.21  and  Fig.22. 


Figure  21.  Velocity  Field  ,  With  V  •  B  =  0  control 


Figure  22.  By  Field.,  With  V  B  =  0  control 


By0=  4  Tesla 


Figure  23.  Comparison  of  Full  MFD  and  Low  Re  Approaches  (red  squares)  for  Normalized  Thrust  and 
Maximum  Velocity  with  Byo=4  T. 

Fig.23  compares  the  Full  MFD  and  Low  Magnetic  Reynolds  Number  Approaches  for  predicting  thrust 
and  maximum  velocity  within  the  accelerator.  At  Byo=4  Tesla  the  two  approaches  agree  very  well. 
However,  note  the  velocity  disturbance  near  x=0  in  the  Full  MFD  results.  It  appears  that  a  shock  wave 
is  propagating  upstream  toward  the  boundary.  Further  inspection  did  show  the  presence  of  a  shock 
wave  moving  upstream.  The  Low  Magnetic  Reynolds  Number  Approach  did  not  and  could  not  show  a 
similar  feature.  The  question  is  -  Should  this  shock  wave  be  physically  present  or  not?  Did  it  start  out 
from  numerical  difficulties  produced  by  the  strong  discontinuity  in  the  imposed  magnetic  field?  The 
calculation  was  repeated  using  as  the  initial  solution  the  converged  solution  from  Low  Magnetic 
Reynolds  Number  Approach.  The  shock  wave  reappeared. 

A  second  test  case  was  simulated  with  a  stronger  magnetic  field.  The  imposed  magnetic  field  was 
Bya  =  11  Tesla  and  the  electric  field  was  Ez0  =  -30,000  V/m.,  again  with  the  same  load  factor  of  2.004. 

All  other  conditions  were  unchanged. 


By0=  1 1  Tesla 


Figure  24.  Comparison  of  Full  MFD  and  Low  Re  Approaches  (red  squares)  for  Normalized  Thrust 
and  Maximum  Velocity  with  By0=ll  T. 

Fig.24  compares  the  Full  MFD  and  Low  Magnetic  Reynolds  Number  Approaches  for  predicting  thrust 
and  maximum  velocity  within  the  accelerator.  Again,  the  presence  of  a  strong  shock  is  observed.  It  has 
moved  to  the  upstream  boundary.  There  are  significant  differences  in  thrust  predicted  by  the  two 
approaches.  At  present,  from  further  testing,  it  is  believed  that  the  Low  Magnetic  Reynolds  Number 
Approach  is  incorrect  and,  at  these  conditions,  a  shock  wave  does  move  upstream  toward  the  entrance 
boundary. 

8.  Inclusion  of  Cesium  Seeding  into  the  Flow  and  the  Park  Model  for  Calculating  Electrical 
Conductivity  Based  upon  Electron  and  Ion  Concentrations  (Item  (9)  in  the  Status  of  Effort  Section 
above). 

Previously  a  numerical  difficulty  was  discussed  in  Sec.2  concerning  the  Park  program  for  calculating 
electrical  conductivity  near  the  wall  within  the  boundary  layer.  This  difficulty  is  now  removed  by 
finding  and  correcting  a  statement  in  the  Park  program.  Simulation  results  were  then  made  for  a  non¬ 
equilibrium  flow. 

An  experimental  Magneto-Hydrodynamic  (MHD)  channel  is  being  tested  at  NASA  Ames  Research 
Center  by  D.W.  Bogdanoff,  C.  Park  and  U.B.  Mehta  to  study  critical  technologies  related  to  MHD 
bypass  scramjet  propulsion.  The  channel  about  a  half  meter  long  contains  a  nozzle  designed  for  Mach 
2  flow,  a  center  channel  section  and  an  accelerator  section.  Magnetic  and  electric  fields  can  be  imposed 
upon  the  flow  within  accelerator  section  for  MHD  acceleration.  The  channel  was  uniformly  2.03cm 
wide. 
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Figure  25.  NASA  MHD  Channel  experiment 
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The  initial  conditions  for  the  nozzle  section  consisted  of  the  flow  at  rest  at  8.0x1 05N/m2  pressure, 
7500°K  temperature  and  an  exit  pressure  low  enough  to  generate  supersonic  flow.  Because  of  flow 
symmetry  only  half  the  channel  was  simulated  in  this  2-D  calculation. 


The  imposed  magnetic  field  was  By—  4  Tesla  and  the  electric  field  was  Ez  — 25,000  V/m.  The  gas  in 

the  accelerator  was  seeded  and  considered  to  be  in  non-equilibrium.  The  following  four  reaction 
seeded  air  chemistry  model  of  Park,  Mehta  and  Bogdanoff  was  used.  The  model  contains  eight  species. 

1)  02  +  M  ^O+O  +  M 

2)  N7  +  O^NO  +  N 

2  B 

3)  0  +  N0<r*N  +  02 

4)  Cs  +  e~  Cs+  +  e~  +  e~ 

B 

The  electrical  conductivity  was  calculated  from  the  species  mole  fractions,  temperature,  density  and 
pressure  of  the  gas,  using  a  program  developed  by  Park.  A  Baldwin-Lomax  turbulence  model  modified 
for  MFD  was  used. 
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Figure  26.  Velocity  profiles,  2-D  accelerator 
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Figure  27.  Maximum  velocity  in  the  axial  direction  within  accelerator  section  for  initial  Cs 
concentrations  of  0.3,  0.6,  1 .2  and  2.4x10 3 


Figure  28.  Normalized  force  in  the  axial  direction  within  accelerator  section  for  initial  Cs 
concentrations  of  0.3,  0.6,  1.2  and  2.4x1  O'3 

Fig.25  shows  the  acceleration  of  velocity  within  the  channel  and  Fig.26  and  Fig.27  show  the  increase 
in  both  axial  velocity  and  force  with  increasing  seeding  concentrations  of  Cs.  Earlier  reported 
difficulties  in  solving  this  problem  have  been  overcome. 

9.  Analysis  of  the  Physics  of  Upstream  Influence  of  the  Flow  about  the  Nose  of  a  Hypersonic  Vehicle 
Concentrations  (Item  (101  in  the  Status  of  Effort  Section  above). 

A  preliminary  study  on  the  upstream  influence  of  the  hypersonic  flow  about  a  sphere  cone  body 
simulating  the  flow  about  the  nose  of  a  hypersonic  vehicle  has  been  started.  We  observed  the 
importance  of  the  upstream  influence  of  flow  governed  by  the  equations  of  magneto-fluid  dynamics 
earlier  within  channels  simulating  the  flow  in  generators  and  accelerators  relevant  to  the  “energy 
bypass  scramjet  engine  concept”.  Now  we  examine  external  hypersonic  flow.  Fig.29  shows  “Mach 
One”  contour  about  a  sphere  cone  body  immersed  in  Mach  1 5  flow.  The  cone  half  angle  is  1 5  degrees. 
The  Mach  One  contour  was  calculated  by  dividing  the  velocity  normal  to  the  body  surface  by  the 
speed  of  sound  and  locating  values  of  unity  .  It  is  essentially  the  location  of  the  bow  shock  wave  about 
the  body,  but  not  exactly  so,  and  separates  the  flow  disturbed  by  the  body  from  the  free  stream. 


Figure  29.  Mach  One  contour  location,  |v„  /c| ,  no  magnetic  field 


Figure  30.  Mach  One  contour  location,  |v„  /c| ,  Interaction  Parameter  Q  =  <JeB%  /  p0u0  =  6 ,  Full  MFD 
Equations 

Fig-30  shows  the  Mach  One  contour  location  with  a  magnetic  dipole  placed  at  the  origin.  Note  the 
increased  standoff  distance.  Fig. 31  shows  the  same  flow  calculated  by  the  Low  Magnetic  Reynolds 
Number  Approach.  The  agreement  of  the  Full  MFD  and  Low  Kri,  MFD  Equations  is  excellent. 


Figure  31.  Mach  One  contour  location,  |v„ /c\ ,  Interaction  Parameter  Q  =  <7eB%  /  p0u0  =  6 ,  Low  Rm 
MFD  Equations 


Fig. 32  shows  the  Mach  One  contour  using  the  “fast  sound  speed”  cf  instead  of  the  usual  speed  of 
sound  c. 
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Magnetic  fields  can  propagate  small  disturbances  with  speed  cf ,  which  can  be  many  times  larger  than 


c .  The  sound  speed  cf  is  directional.  That  given  above  is  for  the  x-direction  of  travel. 


Figure  32.  “Fast  “Mach  One  contour  location, |v„  tcf  | ,  Interaction  Parameter  Q  =  <JeB2ol0  / p0u0  =  6 , 

Full  MFD  Equations 

Note  the  greatly  extended  region  of  disturbed  flow  described  by  the  Full  MFD  Equations.  This 
extended  region  would  not  occur  from  the  Low  Rm  MFD  Equations  because  there  are  not  other  modes 

of  wave  propagation  than  the  usual  sonic  speed. 

An  important  question  is  “what  is  happening  in  the  extended  region?’.  It  appears  that  the  flow  is  slowly 
compressing  and  heating  up  until  it  passes  over  the  bow  shock  wave  where  suddenly  jumps  up  in 
pressure,  density  and  temperature.  Thus,  the  region  just  ahead  of  the  bow  shock  is  not  undisturbed  free 
stream  flow,  but  flow  with  a  higher  temperature  than  that  of  the  free  stream.  Hence,  the  Mach  number 
just  ahead  of  the  bow  shock  wave  is  reduced  from  the  free  stream  Mach  number.  It  is  amazing  that  the 
Low  Rm  MFD  Equations  miss  this  entirely  yet  still  predict  the  bow  shock  location  accurately.  But,  are 

there  other  important  features  of  the  flow  that  the  Low  Rm  MFD  Equations  fail  to  predict?  It  appears 
that  the  heating  of  the  flow  outside  of  the  shock  wave  may  have  a  significant  effect  of  body  surface 
heat  transfer.  This  will  be  explored  subsequently. 


10.  Importance  of  this  Research 

This  research  has  pointed  out  the  simulation  differences  between  the  Full  MFD  and  Low  Reynolds 
Number  Equation  approaches.  The  former  is  a  very  complicated  set  of  equations,  with  non-uniqueness 
issues  and  constraint  satisfaction  difficulties.  However,  it  is  a  more  complete  description  of  the  flow 
about  and  through  bodies  of  Air  Force  interest  than  the  simpler  Low  Reynolds  description,  now  in 
widespread  use.  The  present  status  in  the  development  of  algorithms  for  the  equations  of  magneto¬ 
fluid-dynamics  is  not  yet  mature  enough  to  simulate  flows  about  aerospace  configurations  with 
confidence.  It  is  important  to  continue  this  research  to  avoid  catastrophic  design  errors  caused  by  flow 
simulations  that  do  not  include  all  the  relevant  physics. 
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